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Abstract 

A numerical explicit method to evaluates transient solutions of linear partial dif- 
ferential inhomogeneous equation with constant coefficients is proposed. A general 
form of the scheme for a specific linear inhomogeneous equation is shown. The 
method is applied to the wave equation and the diffuse equation and is investigated 
by simulating simple models. The numerical solutions of the proposed method show 
good agreement to the exact solutions. Comparing with explicit FDM, FDM shows 
the instability by the violation of CFL condition whereas the proposed method is 
always stable irrespective of any time step width. 
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1 Introduction 



Many numerical method have been proposed to solve partial differential equa- 
tions. Many of the proposed method is for spatial modeling. In spite of the 
variety of the spatial modeling, temporal modelings may not have been pro- 
posed so many. 

As for the temporal modeling, when you solve a harmonic oscillating model 
in time, the time differentiation can be replace with iu. When you solve a 
transient model in time, the time differentiation is replaced with finite differ- 
entiation. The finite differentiation can be said to be the most popular method 
for the transient numerical solution. 
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The finite differentiation can be classified roughly into two categories. One is 
the explicit scheme and the rest is the implicit scheme [2] . The explicit scheme 
is easy to program and the amounts of computation for advancing one time 
step is smaller than the implicit method. One of the bad things is the time 
step width must be small enough to satisfy a condition otherwise numerical 
instability may occur. As for the implicit scheme, the computation is generally 
large for one time step but it is always stable for any time step width. From the 
user's point of view the explicit method may be comfortable for programming. 

In the meantime, Green's functions are well used in theoretical physics but 
are not used in the numerical transient problems. This is because if you adopt 
Green's function to solve a certain transient problem you have to integrate all 
over the history in time and space for each time step. Although the Green's 
function may not be suitable for numerical simulation, the Green's function 
is the response function of the basic equation against Dirac's 5 function or 
another functions. This means that the Green's function gives the exact time 
variation of the system against a certain stimulation. 

The reason of the instability of the explicit Finite Difference Method (FDM) is 
caused from the difficulty of the approximation of d/dt with the several order 
of polynomials of At. There is no physical or mathematical inevitability for 
the approximation. Whereas the Green's function is the exact solution of the 
basic equation so the time variation predicted by the Green's function has the 
inevitability. It is a big motivation that the benefit of the Green's function is 
to be integrated into the explicit time advancing to formulate a stable explicit 
scheme irrespective of any time step width. 

The present paper describes a new scheme for transient solutions using Green's 
function to solve inhomogeneous linear differential equation with constant co- 
efficients. Although the scheme is a explicit scheme, it is always stable irre- 
spective of any time step width different from the explicit FDM. 

The preliminary works was applied to electromagnetics and diffuse equations 
[3j[4j[5j. In the present paper the general formulation of the scheme is de- 
scribed. 

In the paper, generalized basic formulation is derived in section [2J Section 
|3] describes specific formulations for wave and diffuse equations. Numerical 
simulation using the proposed method is investigated in section H] to show the 
effectiveness of the presented method. The numerical accuracy and stability 
is discussed in section EJ 



2 



2 Basic Formulation 



Suppose a partial differential equation, 

V t u + V x u = S(x,t), (1) 



where V t ,V x and u are partial differential operator with respect to time t, 
space x = (x, y, z) and a certain physical value, respectively. The partial 
differential operators V t , V x are represented as, 



Qn 



Ql Qra Qn 
l,m,n f 

where a n , bi mn are constant values. The authors are now interested in solving 
Eq.flT]) to obtain the transient solution u for the given source function S{x, t). 

To solve the Eq.flT]) by the proposed method, Fourier transformation is applied 
to the basic equation at first. A physical value can be expressed with integral 
form of Fourier components as, 

oo 

A(M) = 7^372 / MKty^dk, (4) 



where k = (k x , k y , k z ) is the wave number. Substituting Eq.flT]) into Eq.flT]) we 
have, 

V t u(k,t)+JCu(k,t) = S(k,t), (5) 



where 



/C = £ b lmn {ik x )\ik y ) m {ik z )\ (6) 

Z,m,n 



m and S is Fourier components with respect to u and 5*, respectively. 
Suppose a homogeneous equation relating Eq.fl5]), 

(V t + /C)tt(k, t) = C(D)u(k, t) = 0, (7) 
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where D = d/dt. The general solution of Eq.flT]) is obtained as, 



m— 1 
j 1=0 



where m is the multiplicity of the solution Xj, Aj is the solution of C(X) = 
which is the eigen equation of Eq.(J7J). 

Let m be 1, which means there is no multiplicity we have simplified general 
solution of Eq.(jZJ). 



u 



!M) = E<*M (9) 



To obtain explicit scheme for the numerical solution of Eq. (jSJ) we should have 
a recurrence formulation with respect to time step t n . The explicit scheme is 
obtained directly from Eq.® as, 



u 



v k, t n ) — ^2 c j e 



j 

= Cje X: ' t "~ 1 e Xj ^ tn , 
j 

where 

At n = i n - t n _L (10) 

Here the authors introduce a state variable for numerical formulation, 
which is defined as, 

F* (ii) 

Then the authors have explicit scheme for the Eq.flTJ), 



j 

= ^p-i. e ^A. (12) 



General solution of Eq. (jTJ) have been evaluated with a condition of no- multiplicity 
of the eigen values. The derived formulation is recursive with respect to time 
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tn-1 t n 

Fig. 1. Impulse function for the Green's function. 

step. This means u(k, t n ) can be obtained from the state variable of one time 
step before. 

Now the authors consider the general solution of Eq.flT]). The authors introduce 
an impulse function in order to obtain the general numerical solution of the 
inhomogeneous equation (JHJ), 



(t < * n _ x ) 

1 (t n _i < t < t n ) 
(t n < t) 



(13) 



S(k,t) 




t 

Fig. 2. Approximation with the superposition of the Impulse functions. 
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Suppose the Green's function G n {t) of Eq.Q with respect to the impulse 
function which satisfies the equation; 

V t G n {t) + KG n (t) = I(t; t n , At). (14) 

A certain source function S(k, t) can be approximated with superposition of 
the impulse functions as, 

N 

S(k, t)^J2 Sn(K t n - 1/2 )I{t; t n , At), (15) 

n=0 

where the time series S n (k, ^—1/2) is given. 

Equation (j!4p is a linear equation so a superposition of the solutions is also 
its solution. From Eqs. f|T5|) . f|T4l we obtain a superposed solution as, 

N 

u{k, t) « ^(k, t) = J2 S n (K t n )G n (t), (16) 

n=0 

where t > tpj. 

The authors split summation in Eq. ffTB"]) into two parts to obtain an explicit 
scheme for the Eq.([5]), 

N-l 

u N (k, t)=J2 S n (k, t n . 1/2 )G n (t) + S N (k, t N _ 1/2 )G N (t). (17) 

n=0 

If t > tjv, the first term of the right hand of the Eq. (II 71) is the solution of 
Eq.(j7]). So the first term can be described as, 

N-l 

X) ~S n (k,t n „ l/2 )G n (t)=Y;C^\ 

n=0 j 

= Y, F N-i-e Xj{t ~ tN - l} - (18) 
j 

Now we have the explicit scheme to solve Eq.flS]) with the proposed method 
in a general form. 

u N {k,t) =J2 f n-i ■ e x ^- tn -^ + S N (k,t N - 1/2 )G N {t). (19) 
j 
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In the following section we will investigate two concrete examples using the 
proposed scheme. 



3 Formulation Examples 

We present two examples of the proposed formulation for solving partial dif- 
ferential equations which is derived in the present section. First one is a wave 
equation and second is a diffuse equation. They both have constant coeffi- 
cients. 

3.1 Wave Equation 

Here we consider a wave equation of constant coefficients. 



where u is a physical value, x and t are the coordinates in space and time 
respectively, c is the speed of the wave propagation. 5* is a source function and 
is given. We shows the formulation procedure of Eq. (1201) in the following. 

Applying Fourier transform to Eq. ([20]) gives, 



where k 2 = k 2 + k 2 + k 2 . Now we have an ordinary differential equation instead 
of partial differential equation. 

Consider the Green's function of Eq. (l2Tj) against the impulse function (Eq. ffTB"]) ). 
we have an equation, 




(20) 




(21) 




(22) 



Then the eigen equation can be obtained as, 




(23) 



so we have the eigen solutions. 



A = ±ikc. 



(24) 
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According to Eq. (IT9"|) the solution can be described as, 



u N (k, t) = F+_ x ■ e tkc{t - tN - l] + F^_ x ■ e -M*-^-i) 
+S N (k,t N )G N (t), 



(25) 



where F£ and F N are the state variables for A = ike and A = — ike, respec- 
tively. 



Equation ([22]) can be solved analytically as, 



G N (t) = -i 



sin(kcAt N /2) 



k 2 



e ikc(t-t N +At N /2) _ e -ikc(t-t N +At N /2) 



(t > t N ). 



(26) 



Equation (|26|) is the Green's function which satisfies Eq. (l2~2l) . 



We redefine the state variable Fn as, 

F± = • e ±ifccA '- + 5(A;,t iV _ 1/2 )e 



+JkcAt N /2 kcAt N 



(27) 



Then we have the formulation by the proposed scheme. 



u(k, t) 



% 



ikc(t-t N ) rp+ -ikc(t-t N ) tt>- 

e r N — e r N 



(t > t N ). 



(2* 



3.2 Diffuse Equation 



Suppose a diffuse equation with constant diffusion coefficient as, 



^^-cV 2 u(x,t) = S(x,t), 



(29) 



where c is the diffusion coefficient, w(x, t) and S"(x, t) are functions of space x 
and time t. Applying Fourier transformation to the Eq. (I29[) to get the Fourier 
transformed equation. 



du(k,t) 
dt 



+ ck 2 u(k,t) = S(k,t). 



(30) 
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Same as the procedure deriving the formulation for the wave equation we have 
an equation to consider as, 

9G n (t) , _ t2 



+ ck 2 G n (t) = I(t,t n ;At). (31) 



dt 



Equation ( 13T1) can be solved analytically and then we have, 

G N (t) = ^ "y .^^. (32) 



Same as the formulation of the wave equation we have the explicit form of the 
solution as, 

u N (Kt)=F N (k,t N )e~ ckHt ~ tN) , (33) 
where 



F N (k) = F N ' 1 (k)e- CK ' LSlN 

sinhck 2 At N /2 fc2At /2 ~ 



In the next section we will verify the effectiveness of the obtained formulations 
by applying to a simple one dimensional problems. 



4 Numerical Experiments 



A numerical experiments by the proposed scheme are shown here. For reader's 
understanding simple one dimensional examples are shown but the scheme can 
be applied to multi-dimension straightforwardly. 

The source term is given and the same source function is used for both ex- 
amples. The source function adopted here is a Gaussian distribution in space 
and sinusoidal in time. 

f ft < 0) 

S(x,t) = \ (x _ xq)2 ; , (34) 

I e sinw t (t > 0) 



where Xq, Ax are the centre and the half- value width of the Gaussian distri- 
bution, respectively. 
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Table 1 

Common Parameters for the Simulation 



Grid Number 


Length of the system 


Grid Width Ax 


64 


10.0 


10.0/63 



/ \ 
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/ \ 
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(a) Spatial distribution of the source (b) Temporal variation of the source 
function. function. 

Fig. 3. The source function. 
The Fourier transformed Gaussian function S'(k) is analytically evaluated as, 

S(k) = (2nAx 2 ) 1/2 e- k2Ax2/2 e- l]lX0 . (35) 



Figure 3(a) and 3(b) shows the distribution of the function in space and time, 
respectively. The common data in both wave and diffuse equations is shown 
in Table [TJ 

The proposed method is denoted as transient Green method (hereafter TGM) 
for convenience. 



4-1 Wave Equation 



A simple one dimensional numerical simulation is performed in order to in- 
vestigate the characteristics of the TGM. The basic equation to solve is Eq. 
(|21|) and finite difference solution is also calculated for comparison. 

The analytical solution of the basic equation can be obtained and we have, 

,— fc 2 <5^/2— ikxo 



u E (k,t) 



ui - (cky 



— sin ckt — sin coot 
ck 



(36) 



The particular parameters for the wave equation is listed in Table 12 
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Table 2 

Parameters for wave equation. 



Wave Speed (c) 


Time Step Width At 


Frequency (oj/2tt) 


1.0 


0.01 


3.14159 



3 




Fig. 4. Snapshot of Ujv(x, t) and «e(x, i) at t = 1.0. 

The snapshot at i = 1.0 is shown in FigJH Although time advancing is per- 
formed in k space, the figure shows the real space distribution by inverse 
Fourier transformation to help the reader's understanding. Fast Fourier Trans- 
form is used to calculate. The result shows the good agreement with TGM 
and exact solution. 

Figure [5] shows the standard deviation from the exact solutions in k space as 
a function of time step width (At). 



Er(t;At) = /$>(k, *l At) - u E (k, t)) 2 . 



(37) 



In Figj5]the FDM and TGM show error with 0(At 2 ). When the time step 
width At becomes big enough to violate the CFL condition the result of FDM 
shows the numerical instability. While TGM doesn't show any numerical in- 
stability. It seems to be always stable irrespective of any time step width. 
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Fig. 5. Standard deviation of TGM and FDM from the exact solution as a function 
of At. 

4-2 Diffuse Equation 



Same as Sec. I4.1[ numerical simulation is performed with TGM and FDM for 
solving a simple one dimensional diffuse equation problem. The equation to 
solve is a diffuse equation (Eq.( 15Uj) ). The analytical exact solution is obtained 
for this equation, which is: 



u E (k,t) 



e -k 2 5l/2-ikx Q 

(ck 2 ) 2 + ul 

-ck 2 t 



(3* 



oo e " + ck sin u t — u cos cu t . 

The particular parameters adopted in the experiment is listed in Table [3j 
Table 3 

Parameters particularly used for Diffuse Equation. 



Diffuse Coeff. (c) 


Time Step Width At 


Frequency (u/27r) 


3.0 


0.001 


20.0 



The snapshot at t — 0.1 is shown in FigJHJ The solution computed by TGM 
and exact solution indicate good agreement. 



The errors of the TGM and the FDM are shown in FigJTl The error is evaluated 
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_.| Fig. 6. Snapshot of iijv(x, i) and ue(x, t) at t = 0.1. 
form Eq. (1371) same as the case of wave equation. Figure [7] shows numerical 
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Fig. 7. Standard deviation of TGM and FDM from the exact solution as a function 
of At. 

instability is occurred with the violation of CFL condition whereas the results 
of the TGM doesn't show any numerical instability irrespective of any time 
step width. 
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The FDM is the Euler's scheme so the error is O(At). According to the FigfTJ 
TGM has accuracy of 0(At 2 ). 



5 Discussion 



5. 1 Numerical Accuracy 



The numerical error of the proposed scheme potentially contains the discretiza- 
tion errors in time and space. 

The spatial accuracy of TGM is determined by the Fourier transformation. 
If the Fast Fourier Transformation is adopted the data can be transformed 
exactly. This means the transformed data contains no truncation error essen- 
tially. 

The temporal accuracy of the proposed method depend upon the accuracy of 
the integration method of the source term 5(k, t). In the presented scheme 
here (Eqs.(l28l).(l33l)) rectangle mid-point rule is adopted for time integration 
(Eq.( fi~5]) ). So discretization error of Eqs. (l28]) . (l33l is (9 (At 2 ). If you adopt more 
accurate integration scheme the accuracy of the proposed scheme also becomes 
more accurate. 



5.2 Numerical Stability 



The numerical stability of the proposed method is discussed in the section. 

As for the finite difference method the explicit scheme is unstable unless the 
time step width At is small enough to satisfy a condition. The criteria is well 
known as Courant Friedrichs and Lewy (CFL) condition. For instance CFL 
condition of the simple explicit scheme for diffuse equation ( 1291) is; 

2Ar 2 

At < — . (39) 



If the condition is violated the simulation result may be far from the exact 
solution which you want to obtain. The unstableness comes from the difficulty 
in approximating the solution with first or second order polynomials. While 
the proposed method is derived from the superposition of the Green's function 
which are analytical solutions. 
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Now we consider the numerical stability of the TGM solution for Eq.flT]). The 
TGM formulation is given by Eq. (fl2|) . The general solution of Eq.flTJ) is already 
given by Eq.flH]). The state variable F£ also given by Eq.l TTT]) . 

When F 3 n _ x satisfy the equation: 

K_ x = c^~\ (40) 

The Fl can be in TGM formulation as, 

n n—1 ° ) 

= Cje Xjtri . 

So we find the TGM solution is exactly identical with the analytical solution. 
This is always true irrespective of any time step width At. This result indicates 
the TGM formulation is always stable with any time step width as is shown 
in section @] 



6 Conclusion 

A numerical method to compute explicitly for transient solution of linear par- 
tial differential inhomogeneous equation with constant coefficients is proposed. 
It is shown that the proposed method gives the transient solution from the 
state variables at the one time step before explicitly and always numerically 
stable irrespective of any time step width. 
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